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^ , Abstract 
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• An energy functional for orbital based 0{N) calculations is proposed, which 

(N 



depends on a number of non orthogonal, localized orbitals larger than the 
number of occupied states in the system, and on a parameter, the electronic 
chemical potential, determining the number of electrons. We show that the 



^ ■ minimization of the functional with respect to overlapping localized orbitals 

• can be performed so as to attain directly the ground state energy, without 

'x 

■ being trapped at local minima. The present approach overcomes the multiple 

minima problem present within the original formulation of orbital based 0{N) 
methods; it therefore makes it possible to perform 0{N) calculations for an 
arbitrary system, without including any information about the system bond- 
ing properties in the construction of the input wavefunctions. Furthermore, 
while retaining the same computational cost as the original approach, our for- 
mulation allows one to improve the variational estimate of the ground state 
energy, and the energy conservation during a molecular dynamics run. Sev- 
eral numerical examples for surfaces, bulk systems and clusters are presented 
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and discussed. 
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I. INTRODUCTION 



Most electronic structure calculations performed nowadays in condensed matter physics 
are based on a single particle orbital formulation. Within this framework, the ground state 
energy (Eq) of a multi-atomic system is obtained by solving a set of eigenvalue equations. 
Until recently, this has been accomplished by searching directly the eigenstates of the single 
particle Hamiltonian (H), which in general are extended states, e.g. Bloch states in a 
periodic systemB. 

In the last few years, methods for electronic structure (ES) calculations have been intro- 
duced, which are based on a Wannier-like representation of the electronic wave functionsiH^. 
The main motivation for choosing such a representation was the search for methods for which 
the computational effort scales linearly with system size {0{N) methods). Very recently, 
real space Wannier-like formulations were also used to describe the response of an insulator 
to an external electric field!. Within these approaches, a suitably defined total energy func- 
tional (E) is minimized with respect to orbitals constrained to be localized in finite regions 
of real space, called localization regions. The minimization of the energy functional does 
not require the computation of either eigenvalues or eigenstates of H. 

In the absence of localization constraints, one can provei that the absolute minimum of E 
(Eq) coincides with Eq. In the presence of localization constraints, a variational approxima- 
tion to the electronic wave functions is introduced and therefore Eq lies above Eq. However, 
the difference between Eq and Eq can be reduced in a systematic way, by increasing the 
size of the localization regions. We note that localization constraints do not introduce any 
approximation when the resulting localized orbitals can be obtained by a unitary transfor- 
mation of the occupied eigenstates. Therefore the use of localized orbitals is well justified 
for, e.g., periodic insulators, for which exponentially localized Wannier functions can be 
constructed by a unitary transformation of occupied Bloch stateslll. 

The minimization of the functional E with respect to extended states can be easily 
performed so as to lead directly to the ground state energy Eq, without traps at local 
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minima or metastable configuration^. On the contrary, the minimization of E with respect 
to localized orbitals can lead to a variety of minimai'l. In order to attain the minimum 
representing the ground state, information about the bonding properties of the system has 
to be included in the input wavef unctions. This implies a knowledge of the system that may 
be available only in particular cases, and it constitutes the major drawback of the orbital 
based 0{N) method, which has otherwise been shown to be an effective framework for large 
scale quantum simulations^. 

In this paper, we propose a functional for orbital based 0{N) calculations, whose mini- 
mization with respect to localized orbitals leads directly to a physical approximation of the 
ground state, without traps at local minima. This overcomes the multiple minima problem 
present within the original formulationSi and makes it possible to perform 0{N) calculations 
for an arbitrary system, with totally unknown bonding properties. The present formulation 
has also other advantages with respect to the original one. While retaining the same compu- 
tational cost, it allows one to decrease the error in the variational estimate of Eq, for a given 
size of the localization regions, and to improve the energy conservation during a molecular 
dynamics run. 

The novel functional depends on a number of electronic orbitals (M) larger than the 
number of occupied states {N/2) of the A^-electron system, and contains a parameter rj 
determining the total charge. During the functional minimization rj is varied until the 
total charge of the system equals the total number of electrons; thus when convergence is 
achieved, i.e. the ground state is attained, the value of rj coincides with that of the electronic 
chemical potential /i. Once the ground state is obtained for a given ionic configuration, 
the corresponding wave functions and ionic positions can be used as a starting point for 
molecular dynamics simulations, which are then performed at fixed chemical potential. This 
is at variance with conventional ES calculations based on orbital formulations, where is 
always fixed, e.g. by imposing orthonormality constraints. Similar to the present approach, 
0{N) calculations based on a density matrix formulationlii are performed at fixed chemical 
potential. Consistently, the functional describing the total energy does not have multiple 



minima in the subspace of localized density matrices. However, whereas a density matrix 
approach presupposes the use of all the occupied and unoccupied states (i.e. a number 
of states equal to nbasis, where nbasis is the number of basis functions), in our formulation 
only a limited number of unoccupied states needs to be added to the set of occupied states, 
regardless of the basis set size. Therefore the present formulation can be efficiently apphed 
also in computations where the number of basis functions is much larger than the number 
of occupied states in the system (e.g. first principles plane wave calculations). 

The rest of the paper is organized as follows: In section II we present a generalization 
of the original formulation of orbital based 0{N) approaches; we first introduce an energy 
functional which depends on a number of orbitals larger than the number of occupied states, 
and we then discuss its properties and the role of the chemical potential. In section III we 
present the results of tight binding calculations based on the generalized 0{N) method, 
showing that the novel approach overcomes the multiple minima problem, and allows one 
to improve on variational estimates of the ground state properties and on the efficiency of 
molecular dynamics simulations. Conclusions are given in Section IV. 

II. ELECTRONIC STRUCTURE CALCULATIONS AT A GIVEN CHEMICAL 

POTENTIAL 

A. Definition of tiie functional 

We consider the energy functional E defined in Ref. 5, which depends on N/2 occupied 
orbitals, for a N electron system. We generalize E so as to depend on an arbitrary number 
M of orbitals, which can be larger than the number of occupied states N/2. For simplicity, 
we consider a non self-consistent Hamiltonian; however the conclusions of this section are 
easy to extend to self-consistent Hamiltonians. The energy functional is written as: 

M 

E[{0},77,M] = 2 ^ Qij < (t>j\H-v\(t>i > +vN. (1) 

ij=l 
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Here {0} is a set of M overlapping orbitals, H the single particle Hamiltonian, 77 a parameter 
and Q a (M x M) matrix: 

Q = 21 - S. (2) 

S is the overlap matrix: Sij =< > and I is the identity matrix. This definition of the 
Q matrix corresponds to truncate the series expansion of the inverse of the overlap matrix 
to the first order (jV = 1, in the notation of Ref. 5). The charge density is defined as 

M 

p{r) = 2 ^ < (/.j|r >< r\(f)i > Q^j. (3) 

ij=l 

For M — N/2, one recovers the original energy functional for 0{N) calculations. 

We note that the energy functional in Eq. (1) can be expressed in terms of a density 
matrix (3"[{0}] : 

f], M] = 2Tr[{H -r])a] + r]N (4) 
Here the trace is computed over the nbasis functions used for the expansion of the {0}, and 

Before discussing the use of the functional of Eq. (1) within a locahzed orbital formula- 
tion, it is useful to assess some of its general properties. 

(i) E[{0}, 77, M] is invariant under unitary transformations of the type 0^ = S^^i 
where U is a (M x M) unitary matrix. 

(ii) Orbitals with vanishing norms do not give any contribution to the energy functional 
E[{0},77,M]. If the overlap matrix S entering Eq. (1) has (M — M') eigenvalues equal to 
zero, then a unitary transformation U exists, such that {4>'} satisfies the condition: 

< </,^|</,^ >=0, for i = M' + l,...,M. (5) 

Under this condition: 

E[{<f>},rj,M]^E[m,n,M']. (6) 
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We note that if Q is replaced by in the definition of E[{0}, 77, M] (Eq. (|T])), then 
orbitals with a vanishing norm give a non zero contribution to the total energy, since for 

< >^ the eigenvalues of S^^ go to infinity. Therefore the functional E[{0},r7,M], 
with Q replaced by S~^, does not satisfy property (ii). 

(iii) The ground state energy Eq is a stationary point o/E[{0}, 77, M]. In order to prove 
this statement, we consider the following set of orbitals {4>^}: 

10° > = \x^> for t = l,N/2 

|0 > for i = N/2 + l,M (7) 

where \xk > are the nbasis eigenvectors of H with eigenvalue e^. Hereafter we assume that 

< XklXk >= 1 and efc < e^+i. The set {0°} fulfills Eq.(|), and therefore E[{0°},r/,M] = 
Ei[{(f)^},ri, N/2] = Eq. In addition, the set {(jP} is a stationary point of E[{0},r7,M], since 
5E/(50fc|{0O} = 0, where 

5E 

— = 4^[(i/-r^)|0, > (g,fc) - 10, X 0,|(i/-r^)|0, >]. (8) 

(iv) The stationary point Eq is a minimum ci/E[{0}, 77, M] if rj is equal to the electronic 
chemical potential /i. We will only consider electrons at zero temperature, and therefore we 
choose /z such that eAr/2 < ^ < eAr/2+i- This property will be proved in the next section. 

B. Role of the chemical potential 

Before giving a proof of property (iv) stated in section II. A, we discuss a simple example 
which is useful to illustrate the role played by rj in the minimization of the energy functional 
E. For this purpose, we evaluate the functional E[{0}, r], M] for a set of M eigenstates of the 
Hamiltonian. In particular, we choose a set {0} such that |0j >= ai\xi >, with arbitrary 
Qi. In this case the energy functional becomes: 

M 

E[{a}, r^, M] = 2 ^(e, - r/)(2 - a^,)ai + r/iV (9) 

i=l 
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As illustrated in Fig. 1, the function (e^ — 77) (2 — af)a'f has a minimum at Oj = if > 77, 
and a minimum at = 1 if < r/. Thus the functional E[{a}, t], M] has a minimum for a 
set {a°} such that a° = 1 if ej < r], and = if ej > rj. At the minimum, Eq. becomes 

M' 

E^in = 2^e, + r/(Ar-2M'), (10) 

i=l 

where eM' < 1] < eu'+i and the total charge of the system is ptot = 2X]i^i(2 — = 2M'. 
We can now choose rj so that p is equal to the actual number of electrons in the system. 
This is accomplished by setting eN/2 < rj < eN/2+1, i-e. by choosing t] equal to the electronic 
chemical potential p. We then have ptot = 2M' = and and -Emin = Eq. 

In order to give a general proof of property (iv) (section II. A), we show that the Hessian 
matrix (h) of the functional E[{0},?7, M] at the ground state is positive definite, ii r] = p. 
The computation of the eigenvalues of h follows closely the procedure used in Ref. ^ to 
calculate the Hessian matrix of E[{0}, 77, A^/2] at the ground state. Since the functional 
E[{0}, 77, M] is invariant under unitary rotations of the {0}, we can write a generic variation 
of the wave function with respect to the ground state as 

10° > =\Xi> +1 Ai > for z = 1, N/2 

|0>+|A, > for i = A^/2 + l,M (11) 

where 

"^basis 

|A. >= E 4x1 > . (12) 

1=1 

By inserting Eq. ( pT]) into Eq. (1), it is straightforward to show that the first order term 
in the {c} coefficients vanishes for any value of the parameter 77, consistently with property 
(iii) stated in section II. A. The remaining second order term can be written as follows: 

E^'^ = J: E 2[e„-6.](4)^+E8[r/-^^^][^(c;. + c^-)f + 

i=lm=N/2+l ij=l 

M "basis 

E E 4[e^-7^](0^ (13) 

i=N/2+l m=N/2+l 

8 



The eigenvalues 2[em — u] are independent of rj and always positive, whereas the eigenvalues 
8[ri — (ej + ej)/2] and 4(6^ — rj] are positive, if and only if rj coincides with the chemical 
potential fi. This proves property (iv) of section II. A. 

III. 0(iV) CALCULATIONS WITH OVERLAPPING LOCALIZED ORBITALS 
A. Localization of orbitals and practical implementation 

We now turn to the discussion of the functional defined in section II. A within a localized 
orbital formulation. The use of localized orbitals is a key feature to achieve linear system- 
size scaling! calculations. Orbitals are constrained to be localized in appropriate regions of 
space, called localization regions, i.e. they have non zero components only inside a given 
localization region, whereas they are zero outside the localization region. The choice of the 
number of localization regions and of their centers is arbitrary. In the calculations that will 
be discussed in the next sections, we chose a number of localization regions equal to the 
number of atoms, each centered at an atomic site (/) . We then associated an equal number 
of localized orbitals (n^) to a localization region, e.g. two and three localized orbitals for 
M = N/2 and M = 3N/2, respectively. 

We will present electronic structure calculations and molecular dynamics simulations of 
various carbon systems, carried out within a tight binding approach. We adopted the TB 
Hamiltonian proposed by Xu et al.0'0, which includes non zero hopping terms only between 
the first nearest neighbors. In a tight-binding picture, a localization region centered on the 
atomic site / can be identified with the set {LRj} of atoms belonging to the localization 
region. Atoms are included in {LRj}, if they belong to the nearest neighbor of the center 
atom. Then, the localized orbital |0f >, whose center is the Jth atom, is expressed as 



J€{LBi} I 

where \aji >'s are the atomic basis functions of the atom J and the index / indicates the 
atomic components {s,Px,Py or p^). In our computations, the generalized energy functional 




L 




(14) 
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was minimized with respect to the locahzed orbitals {0^} by performing a conjugate gradient 
(CG) procedure, both for structural optimizations and molecular dynamics simulations. For 
some calculations it was necessary to use a non zero Hubbard-like term^l to prevent unphys- 
ical charge transfers. In this case the line minimization required in a CG procedure reduces 
to the minimization of a polynomial of eighth degree in the variation of the wavefunction 
along the conjugate direction. We performed an exact line minimization by evaluating the 
coefficients of the polynomial, and by solving iteratively for the polynomial roots. 

B. The multiple minima problem 

As mentioned in the introduction, the major drawback of the original formulation of 
orbital based 0{N) calculations is the so called multiple minima problem. Experience has 
shown that the minimization of E[{(j)},ri, N/2] with respect to localized orbitals usually 
leads to a variety of minimal0, and that the physical properties of the minimum reached 
during a functional minimization depend upon the choice for the input wave functions. If 
the input wave functions are constructed by taking advantage of bonding information about 
the ground state, then a minimum representing a physical approximation to the ground state 
may be reached, after an iterative minimization. On the contrary, if no information on the 
ground state is included in the localized orbitals from the start, the functional minimization 
usually leads to a local minimum, which is characterized by an unphysical charge density 
distribution. 

This is illustrated for a particular case in Table I and Fig. 2, where we present the results 
of a series of tight binding (TB) calculations using localized orbitals, for a 256 carbon atom 
slab. The slab, consisting of 16 layers, represents bulk diamond terminated by a C(lll)-2 
X 1 Pandey reconstructed surface on each side. We considered localization regions (LRs) 
extending up to second neighbors (N/i=2). We performed conjugate gradient minimizations 
of the electronic structure using two localized orbitals per LR (ns=2), which correspond to 
the case M = N/2 in Eq. (1), i.e. to the original formulation of 0{N) calculations. These 
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minimizations were carried out by starting from different wave function inputs. Tlie only 
calculation which lead to a physical minimum was the one started with orbitals containing 
symmetry information about the system, as shown by comparing the results of Fig. 2C with 
those of direct diagonalization, reported in Fig. 3B. The other calculations lead to unphysical 
minima: when starting with a totally random input (Figs. 2A), we found a local minimum 
with charged sites, located predominantly in the surface layers and in the middle of the 
slab. When starting from an atom by atom input (Fig. 2B) we obtained a local minimum 
corresponding to two differently charged surfaces, one positively and the other negatively 
charged. 

The local minima problem present in the original 0{N) formulation can be illustrated 
with a simple one dimensional modehlll We consider a linear chain with N^ite sites and 2 
iVgite electrons in a uniform electric field of magnitude F, with Hamiltonian: 

H = [Egap\eK >< exl - FK{\eK >< ckI + \9k >< 9k\)]- (15) 

K=l 

Here |ex > and \gK > are the highest and the lowest level of the isolated site K, respectively, 
and Egap is the splitting between these two levels. Since the hopping terms between different 
sites are set at zero, {ex > and {qk > are also eigenf unctions of the linear chain Hamiltonian. 
We now study the ground state of the system as a function of the electric field F. If 
< F < Egap/ (A^sitc — 1)5 the total energy of the system is minimized by the set of orbitals 
10° > given by: 

|0°>=|(7. > for z = l,iV3ite. (16) 

If Egap/ {Nsitc — 1) < F < Egap/ {Nsitc — 2), the eigenvalue of \gi > is higher than that of 
kAfsite ^^'^ therefore the total energy of the system is minimized by the following set of 
orbitals >: 

10° > = \gi+i > for i = 1, N^ite - 1, 

\ei> ioT i = Nsite- (17) 
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In both cases, the total energy of the hnear chain system can be obtained exactly within a 
localized orbital picture, by considering iVsite LRs centered on atomic sites, which extend up 
to the first neighbors of a given site. 

We first describe the total energy of the system with the functional E [{</)}, 77, A^/2]. 
Within this framework, the set |0° > which minimizes E[{0}, 77, in the presence of 
a small field, i.e. when < F < Egap/ {Nsitc — 1), is also a local minimum of E[{0}, 77, A^/2] 
in the presence of a large field, i.e. when Egap/{Nsitc — 1) < F < Egap/{Nsitc — 2). This can 
be easily seen from the second order expansion of E[{(l)},ri,N/2\ around the set of 

orbitals defined in Eq. (|I6|) : 

E^'^=i: E nE,ap-F{m-^)]{elf+Y: E ^Iv + F'-^][^{g] + gj)]' , (18) 

i=l m£{LRi} i=l je{LRi} ^ 

where g]^ and e}^ are the projections of the vector |0j > — 10° > on the state \gk > and \eK >, 
respectively. If the orbital are extended, the difference (m—i) can be as large as (Ai'sitc— 1) and 
the eigenvalues [Egap-F{m-i)] can be negative when Egap/lA^site-l) < F < Egap/{Nsitc-'2). 
However, if the orbital are localized the difference {m — i) is smaller than (A^site ~ 1) and 
the eigenvalues [Egap — F{m — i)] are always positive also for Egap/{Nsite — I) < F < 

Egap/iNsite — 2). 

We now turn to a description of the total energy of the linear chain system with the 
functional Fi[{(f)}, fi, M], where M is larger than the number of occupied states N/2, e.g. 
M = 2iVsite. It is straightforward to show that contrary to a description with E[{0}, rj, N/2], 
when using E[{0}, /i, M] the set of orbitals of Eq. ( |16|) is not a local minimum of the system 
in the presence of a large field. Indeed, according to Eq. ([TsD , the second order expansion 

^(2) 

is now given by: 

-^sitc -^site A _|_ 1 

E''^=T. E ^[E,ap-F{m-^)]{elr+Y. E + F ^][-={g] + g^f + 

i=l m£{LR,} i=l jG{LRi} ^ V ^ 

E E 4[E,-Fm-^](^j^ (19) 
Here the LOs with indices i and i + A^^site are assigned to the localization region {LRi}. 
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Both within an extended and a locahzed orbital picture, the eigenvalue 4[£'p — FNgHe — A*] 
is negative when Egap/iN^ite - 1) < F < Egap/iN^ite - 2). 

This simple model shows that the extremum properties of the functional E[{0}, 77, N/2] 
and E[{0},//, M] are in general different, and in particular that local minima of 
E[{0}, r], are not necessarily so for E[{0},/x, M]. This suggests that the use of the 
functional E[{0},/x, M] can overcome the multiple minima problem encountered within a 
formulation based on E[{0}, 77, A^/2]. This simple model suggests also the reason why the 
multiple minima problem should be overcome: the presence of the global variable /i, to- 
gether with the augmented variational freedom of extra orbitals added to the definition of 
the functional, can account for global changes taking place in the system. 

C. Overcome of the multiple minima problem 

We now present a series of numerical examples, showing that the minimization of the 
generahzed functional E[{0},77,M] (Eq. (1)) with respect to localized orbitals can be per- 
formed without traps at local minima, as indicated by the simple model discussed in the 
previous section. We performed calculations for various carbon systems (bulk solids, sur- 
faces, clusters and liquids), by using again LRs extending up to second neighbors (Nfi—2). 
We considered three LOs per site (n^^S), i.e. M — 3N/2 in Eq. (1). In all cases, using ng—S 
was sufficient to overcome the multiple minima problem present in the original formulation. 
We note that the generahzed functional, although it includes a number of localized orbitals 
larger than the number of occupied states, still allows one to carry out electronic minimiza- 
tions and molecular dynamics simulations with a computational effort scaling linearly with 
system size. 

In Fig. 4, we show the energy and the charge per atom during a conjugate gradient 
minimization of E[{0},?7,M], for a 256 carbon atom slab, starting from a totally random 
input. The system is the same as the one studied in the previous section with ns=2. The 
minimization was started with r) — 20 eV; the parameter was then decreased every 20 
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iterations, and finally set at 3.1 eV, which corresponds to the value of the chemical potential. 
As discussed in section II. B, for a given 77 the integral of the charge density converges to a 
value which corresponds to filling all the orbitals with energies smaller than r]. For example, 
for 77 = 20 eV the total charge per atom is equal to 6, i.e. all the 3N/2 orbitals are filled. 
Eventually, when rj = /j, the total charge becomes equal to the number of electrons in the 
system. The way 77 is varied during a minimization is not unique; however the final value 
of 77 must be always adjusted so as to obtain the correct charge in the system. It is seen 
in Table I that all the minimizations with Ug—S converge to the same value, irrespective of 
the input chosen for the wave function. This value corresponds to a physical minimum, as 
shown in Fig. 3 where we compare the charge density distribution with that obtained by 
direct diagonalization. 



D. Improvement on variational estimates of the ground state properties 

The use of the generalized functional and LOs not only overcomes the problem of multiple 
minima, but it also improves the variational estimate of £'0, for a given size of the LRs. This 
is shown is Table II and III, where we compare the results of calculations using the same LRs 
but different number of orbitals (ns=2 and 3), for various carbon systems. The improvement 
is particularly impressive in the case of Ceo, where we also performed an optimization of the 
ionic structure. The error on the cohesive energy is decreased from 3 to 1.5 % by increasing 
Us from 2 to 3. Most importantly the optimized ionic structure obtained with ns=3 is in 
excellent agreement with that obtained with an extended orbital calculation. We note that 
localization constraints introduce a symmetry breaking in the system, i.e. LOs do not satisfy 
all the symmetry properties of the Hamiltonian eigenstates. In Ceo the symmetry breaking 
is large when using ns—2] the deviation of the double and single bond lengths with respect 
to their average values are 3.5 and 6.3 %, respectively. On the contrary, in the optimized 
geometry obtained with ns=3 the symmetry breaking is very small (0.1 and 0.5 %, for the 
double and single bonds, respectively), compared to the icosahedral structure. 
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When using ns=2, the ground state LOs are nearly orthonormali, whereas minimizations 
with ns=3 yield overlapping LOs. Indeed when using ns=3, at the minimum the overlap 
matrix S has 2ns eigenvalues close to 1 and n^ eigenvalues close to 0, and this condition can 
be satisfied with a non diagonal S matrix. We define a quantity measuring the orthogonality 
of the orbitals as = iEfj=Mj - Sijf)/M. In the case of Cgo, is 2.5 lO^^^ and 0.17 
for Yis=2 and ns=3, respectively. We also note that for various systems, the centers of the 
LOs < >, defined as 

>- <0^|0^> ' ^^^^ 

were always found to be located at distances shorter than one bond length from the center 
of their own LRs, when using ns=3. In the case of tls=2, we instead found cases, e.g. the 
Geo molecule, where some orbitals were centered far from their atomic sites and close to the 
border of their LRs. 

E. Molecular dynamics simulations 

In order to investigate the performances of the generalized functional (Eq. 1) for molec- 
ular dynamics (MD) simulations, we carried out MD runs for liquid carbon at low density 
(2 gr cm~^) and at 5000 K. We used a 64 atom cell with simple cubic periodic boundary 
conditions and only the F point to sample the BZ. We used a cutoff radius of 2.45 A for the 
hopping parameters entering the TB Hamiltonian and for the two body repulsive potentiallil 
and U = 8 eV. In the case of 1-C it was necessary to add an Hubbard like term to the Hamil- 
tonian, in order to prevent unphysical charge transfers during the simulations. Equilibration 
of the system was performed in the canonical ensemble by using a Nose thermostat^. 

Within the original 0{N) approach, MD runs for 1-C were found to be particularly 
demanding from the computational point of view, since they required many iterations (Niter) 
per ionic move (e.g. Niter=300 for At=30 a.u.), in order to minimize the energy functional^. 
Most importantly, during the simulation the system could be trapped at a local minimum, 
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evolve adiabatically from that minimum for some time, and suddenly jump to another 
minimum lower in energy. This shows up as a spike in the constant of motion of the system 
(Econst), as can be seen in the line (c) of Fig. 5, which displays Econst for a run performed with 
ns=2. Because of local minima, a perfect conservation of energy could never be achieved 
with ns—2, even by increasing Niter to a very large number. 

When MD runs are performed with n5=3, the problem of local minima is overcome; 
furthermore a significant improvement in the conservation of energy can be achieved at the 
same computational cost as simulations with ns=2. This is seen in Fig. 5 by comparing lines 
(b) and (c). When the generalized functional is used, the accuracy of the energy conservation 
during a MD run is related only to the convergence of the electronic minimization scheme: 
a good conservation of energy can be obtained just by increasing Niter- This is shown by the 
line (a) in Fig. 5. We note that the behavior of Econst observed for all the simulations was 
not affected by the presence of the thermostat. This was checked by repeating all MD runs 
with three different masses (Q^) for the Nose thermostat (Qs=l,4,100 in the same units). 
The structural properties of 1-C computed from the MD runs with ns=3 showed a very good 
agreement with those previously obtained with ns=2. 

IV. CONCLUSIONS 

We have presented a generalization of orbital based 0{N) approaches, which relies upon 
a novel functional, depending on a number of localized states larger than the number of 
occupied states, and on a parameter which determines the total number of electrons in the 
system. We have shown that the minimization of this functional with respect to localized 
orbitals can be carried out without traps at local minima, irrespective of the input chosen 
for the wave functions. In this way, the multiple minima problem present in the original 
formulation is overcome, and 0{N) computations can be performed for an arbitrary system, 
without knowing any bonding properties of the system for the calculation input. We have 
also presented a series of tight binding calculations for various carbon systems, showing that 
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the generalized 0{N) approach allows one to decrease the error in the variational estimate 
of the ground state properties, and to improve energy conservation, i.e. efficiency, during 
a molecular dynamics run. This can be accomplished at the same computational cost as 
within the original formulation. At variance from 0{N) density matrix approaches, our 
formulation requires that only a limited number of unoccupied states be included in the 
energy functional, regardless of the basis set size. Therefore the present formulation can 
be efficiently applied also in computations where the number of basis functions is much 
larger than the number of occupied states in the system (e.g. first principles plane wave 
calculations) . 
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FIGURES 



Fig. 1 Plot of the function f{ai,rj) — (ej — 77)0^(2 — a?) for a positive and a negative 
value of (ej — 77) . 

Fig. 2 Differential atomic charge (Ap) on each atomic site of a 256 carbon atom slab. 
The slab, consisting of 16 layers, represents bulk diamond terminated by a C(lll)-2 x 1 
Pandey reconstructed surface on each side. The ionic index indicates individual atomic sites 
belonging to the slab, which are ordered layer by layer, starting from the uppermost surface. 
The arrow indicates the slab center. Apx = Pk — where px = '2J2ij=iJ2i < (pilc^Ki > 
Qij < 0£Ki\(l>j >, p^ — 4:, and K is the atomic site. In panels A, B and C we show the results 
of calculations performed with two orbitals per atomic site, and with the three different 
wave function inputs listed in Table I, respectively. 

Fig. 3 Differential atomic charge (Ap) on each atomic site for the same system as in 
Fig. 2. The ionic index is the same as in Fig. 2. In the upper panel we report the results of 
a calculation carried out with three orbitals (n^) per atomic site, and with a totally random 
input for the initial wave functions (see Table I). Contrary to the calculation started from 
a totally random input and performed with ns=2 (see Fig. 2A), the calculation with Us—S 
gives a ground state charge density very close to that obtained by direct diagonalization, 
shown in the lower panel. 
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Fig. 4 Total energy, Etot, (upper panel) and total charge (lower panel) per atom, as a 
function of the number of iterations, for an electronic minimization of the same system as in 
Figs. 1,2 and Table 1. Etot= r], M] (see text). The minimization was carried out with 

three states per atom (n^^S) and was started from a totally random input. The chemical 
potential (ry) was varied from 20 to 3.1 eV during the minimization. The final value of rj 
was chosen so that the total charge eventually be equal to the number of electrons in the 
system. In the upper panel, the inset shows Etot as a function of 500 iterations, converging 
to the value reported in Table I, and indicated as a dotted line. 

Fig. 5 Energy per atom (Econst) as a function of the simulation time (t) for constant tem- 
perature (T) molecular dynamics (MD) simulations of liquid C. Econst=Ekin + E[{0},?7,M] 
+ Ethrms, where Ekin is the ionic kinetic energy, E[{0}, 77, M] is the ground state value of 
the electronic energy functional (see text) and Ethrms is the sum of the potential and kinetic 
energies associated to the Nose' thermostat. The LRs extend up to second neighbors (N/j=2, 
amounting on average to 18 atoms per LR). Lines (a) and (b) correspond to MD runs with 
three states per atom (n^^S), whereas line (c) corresponds to a simulation with ns—2. The 
time step used in the three MD runs was 30 a.u.(0.73 fs). At each step, the electronic struc- 
ture was minimized by a conjugate gradient procedure with a fixed number of iterations 
(Niter)- The simulations represented by lines (b) and (c) require the same computational 
cost. 
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TABLES 



Wave Function input 


Ec K=2] 


Ec K=3] 


Totally random 


6.837 


6.978 


Atom by atom 


6.721 


6.978 


Layer by layer 


6.930 


6.978 



TABLE I. Cohesive energy Ec (cV) of a 256 carbon atom slab. The slab, consisting of 16 layers, 
represents bulk diamond terminated by a C(lll)-2 x 1 Pandey reconstructed surface on each side. 
Ec was obtained by performing localized orbital calculations with two and three states (n^) per 
atom (see text), and with three different inputs for the starting wave functions. Totally random 
input: the wave function expansion coefficients (C}^, see Eq. (14)) on each site of a localization 
region (LR) are random numbers, and orbitals belonging to the same LR are orthonormalized at 
the beginning of the calculations. Atom by atom input: each orbital has a non zero Cji only on 
the atomic site to which it is associated, and for each atomic site this coefficient is chosen to be 
the same. Layer by layer input: each orbital has a non zero Cji only on the atomic site to which 
it is associated, and the value of this coefficient is chosen to be the same for each equivalent atom 
in a layer. In the case of atom by atom and layer by layer inputs, the initial wave functions are 
an orthonormal set. The calculations were performed with 77=7.5 eV and ?7=3.1 eV for ns=2 and 
ns=3, respectively, and with LRs extending up to second neighbors (Nft=2, amounting at most to 
17 atoms per LR). The value for Ec obtained by direct diagonalization is 7.04 eV. (See also Fig. 1). 
The highest occupied and lowest unoccupied eigenvalues are 2.85 and 3.42 eV, respectively. In all 
calculations the Hubbard like term was set at zero. 
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Physical properties 


Cohesive Energy/atom 


Double-bond distance 


Single-bond distance 


L0[Nh=2, n,=2] 


6.69 (6.89) 


1.358-1.407 


1.420-1.512 


L0[Nh=2, n,=3] 


6.81 (6.91) 


1.386-1.388 


1.445-1.453 


Extended Orbitals 


6.91 


1.393 


1.440 



TABLE II. Cohesive energy (eV) and length (A ) of the double and single bonds of Cqo, as 
obtained from structural optimizations using localized (LO) and extended orbitals. In all calcula- 
tions the Hubbard like term was set at zero. For comparison, cohesive energies obtained by direct 
diagonalization are given in parentheses. Computations with LO were performed by including 
two shells in a localization region (N/j=2, amounting to 10 atoms per localization region), and by 
considering two and three orbitals (n^) per atom (see text). 
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Crystal structure 


Diamond (ro = 1.54 A) 


2D-Graphite (ro = 1.42 A) 


ID-Chain (ro = 1.25 A) 


Ee [N,,=2, n,=2] 


7.16 


7.09 


5.62 


Ec [N;,=2, n,=3] 


7.19 


7.12 


5.67 


Ec [Nft=oo] 


7.26 


7.28 


5.93 



TABLE III. Cohesive energy Eg (eV) of different forms of sohd carbon computed at a given 
bond length tq. The calculations were performed with supercells containing 216, 128 and 100 
atoms for diamond, two-dimensional graphite and the linear chain, respectively. In calculations 
with localized orbitals we used 2 and 3 orbitals per atom (n^^, see text). The LRs included two 
shells of neighbors (N/i=2), amounting to 17, 10 and 5 atoms per LR in the case of diamond, 
two-dimensional graphite and the linear chain, respectively. 
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